##################################################################
# Nadiya Kostyuk
# Sensitivity Analysis
# R version: 4.1
# 11/19/25
##################################################################

rm(list=ls())

## Install & load packages (all at once)
list.of.packages <- c('dplyr',  'ggplot2', 'broom')
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]; if(length(new.packages)){install.packages(new.packages,dependencies=TRUE)}
lapply(list.of.packages, require, character.only = TRUE)


# setting directory: 
setwd("")


# loading data
load('data/data_final.RData')

################################
# Maneuver versus Attrition
# Hypotheses
##################################

# subsetting dataset: 
data_d_contested_80 = data_final %>% 
  filter(territorial_control_80_1dlag == 'Contested')


# main model:
model <- glmer(
  disruption_intentional_3indicators ~
    disruption_intentional_3indicators_1dlag +
    as.factor(war_phases_1dlag) +
    log_rus_mil_t_1dlag +
    holiday_dummy +
    as.factor(WEEKEND) +
    (1 | ADM1_NAME),
  data = data_d_contested_80,
  family = binomial()
)
vars_to_track <- grep("war_phases_1dlag", names(fixef(model)), value = TRUE)
vars_to_track


# Function to run model and extract results for selected variables
run_main_model_extract <- function(df) {
  model <- glmer(
    disruption_intentional_3indicators ~
      disruption_intentional_3indicators_1dlag +
      as.factor(war_phases_1dlag) +
      log_rus_mil_t_1dlag +
      holiday_dummy +
      as.factor(WEEKEND) +
      (1|ADM1_NAME),
    data = df,
    family = binomial()
  )
  
  # Extract results for selected variables
  tidy(model) %>%
    filter(term %in% vars_to_track) %>%
    dplyr::select(term, estimate, std.error, p.value)
}

# Sensitivity function

run_sensitivity_sim <- function(df, flip_fraction) {
  df_sim <- df
  
  # Identify unintentional outages
  unint_idx <- which(df_sim$disruption_unintentional_3indicators == 1)
  # flip only those unintentiional outages to intentaionl
  n_flip <- round(length(unint_idx) * flip_fraction)
  flip_idx <- sample(unint_idx, n_flip)
  df_sim$disruption_intentional_3indicators[flip_idx] <- 1
  
  # Run model and extract results
  model_results <- run_main_model_extract(df_sim)
  model_results$flip_fraction <- flip_fraction
  return(model_results)
}
rep_sensitivity_sim = function(df, flip_fraction, n){
  lapply(seq_len(n), function(i){run_sensitivity_sim(df, flip_fraction)})
}


set.seed(123456)
flip_rates <- seq(0, 0.5, by = 0.05) 
results_list <- lapply(flip_rates, function(f) 
  rep_sensitivity_sim(data_d_contested_80, f, 100)) 
results_df <- bind_rows(results_list)
results_df_sum = results_df %>%
  filter(flip_fraction > 0) %>%
  group_by(term, flip_fraction) %>%
  mutate(
    estimate = mean(estimate), 
    std.error = sqrt(var(estimate) + mean(std.error^2))
  ) %>%
  ungroup()
results_df_sum0 = results_df %>%
  filter(flip_fraction == 0) %>%
  group_by(term) %>%
  slice(1)
results_df_sum_final = bind_rows(results_df_sum0, results_df_sum)
sort(unique(results_df_sum_final$term))
# saving results from sensitivy analysis: 
#save(results_df_sum_final, file = 'analysis/output/sensitivity_sim_data.RData')


#####################################
# loading data: 
load('analysis/output/sensitivity_sim_data.RData')
# renaming variables for the plot: 
results_df_sum_final = results_df_sum_final %>% 
  mutate(term = case_when(
    term == 'as.factor(war_phases_1dlag)Rus Maneuver Defensive' ~ 'Rus Defensive Maneuver',
    term == 'as.factor(war_phases_1dlag)Rus Maneuver Offensive' ~ 'Rus Offensive Maneuver'
  )
  
  )

ggplot(results_df_sum_final, aes(x = flip_fraction, y = estimate, color = term)) +
  geom_line(size = 1) +
  geom_ribbon(aes(ymin = estimate - 1.96 * std.error, ymax = estimate + 1.96 * std.error, fill = term),
              alpha = 0.2, color = NA) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(
    title = "Sensitivity of Coefficients to Reclassification of Outage Intentionality",
    x = "Fraction of Unintentional Outages Reclassified as Intentional",
    y = "Coefficient Estimate"
  ) +
  theme_minimal() +
  theme(legend.title = element_blank())
#ggsave(file = 'analysis/output/figures/sensitivity_sim_plot.pdf', width = 10, height = 6)


################################
# Punishment Hypothesis
##################################

# subsetting dataset: 
data_d_no_ruscontrol_80 = data_final %>% 
  filter(territorial_control_80_1dlag != 'Russian occupation')
dim(data_d_no_ruscontrol_80)


# main model:
model <- glmer(
  form = disruption_intentional_3indicators  ~
    disruption_intentional_3indicators_1dlag
  + as.factor(territorial_control_80_1dlag)
  + log_rus_mil_t_1dlag
  + holiday_dummy
  + as.factor(WEEKEND)
  + (1|ADM1_NAME), 
  data = data_d_no_ruscontrol_80, 
  family = binomial()
)
summary(model)
vars_to_track <- grep("territorial_control_80_1dlag", names(fixef(model)), value = TRUE)
vars_to_track


# Function to run model and extract results for selected variables
run_main_model_extract <- function(df) {
  model <- glmer(
    disruption_intentional_3indicators ~
      disruption_intentional_3indicators_1dlag +
      as.factor(territorial_control_80_1dlag) +
      log_rus_mil_t_1dlag +
      holiday_dummy +
      as.factor(WEEKEND) +
      (1|ADM1_NAME),
    data = df,
    family = binomial()
  )
  
  # Extract results for selected variables
  tidy(model) %>%
    filter(term %in% vars_to_track) %>%
    dplyr::select(term, estimate, std.error, p.value)
}

# Sensitivity function

run_sensitivity_sim <- function(df, flip_fraction) {
  df_sim <- df
  
  # Identify unintentional outages
  unint_idx <- which(df_sim$disruption_unintentional_3indicators == 1)
  # flip only those unintentiional outages to intentaionl
  n_flip <- round(length(unint_idx) * flip_fraction)
  flip_idx <- sample(unint_idx, n_flip)
  df_sim$disruption_intentional_3indicators[flip_idx] <- 1
  
  # Run model and extract results
  model_results <- run_main_model_extract(df_sim)
  model_results$flip_fraction <- flip_fraction
  return(model_results)
}
rep_sensitivity_sim = function(df, flip_fraction, n){
  lapply(seq_len(n), function(i){run_sensitivity_sim(df, flip_fraction)})
}


set.seed(123456)
flip_rates <- seq(0, 0.5, by = 0.05) 
results_list <- lapply(flip_rates, function(f) 
  rep_sensitivity_sim(data_d_no_ruscontrol_80, f, 100)) 
results_df <- bind_rows(results_list)
results_df_sum = results_df %>%
  filter(flip_fraction > 0) %>%
  group_by(term, flip_fraction) %>%
  mutate(
    estimate = mean(estimate), 
    std.error = sqrt(var(estimate) + mean(std.error^2))
  ) %>%
  ungroup()
results_df_sum0 = results_df %>%
  filter(flip_fraction == 0) %>%
  group_by(term) %>%
  slice(1)
results_df_sum_final = bind_rows(results_df_sum0, results_df_sum)
sort(unique(results_df_sum_final$term))
# saving results from sensitivy analysis: 
#save(results_df_sum_final, file = 'analysis/output/sensitivity_sim_data_territories.RData')


#####################################
# loading data: 
load('analysis/output/sensitivity_sim_data_territories.RData')
# renaming variables for the plot: 
results_df_sum_final = results_df_sum_final %>% 
  mutate(term = case_when(
    term == 'as.factor(territorial_control_80_1dlag)Ukrainian control' ~ 'Punishment in Ukr Territories'
  )
  
  )

ggplot(results_df_sum_final, aes(x = flip_fraction, y = estimate, color = term)) +
  geom_line(size = 1) +
  geom_ribbon(aes(ymin = estimate - 1.96 * std.error, ymax = estimate + 1.96 * std.error, fill = term),
              alpha = 0.2, color = NA) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(
    title = "Sensitivity of Coefficients to Reclassification of Outage Intentionality",
    x = "Fraction of Unintentional Outages Reclassified as Intentional",
    y = "Coefficient Estimate"
  ) +
  theme_minimal() +
  theme(legend.title = element_blank())
#ggsave(file = 'analysis/output/figures/sensitivity_sim_plot_territories.pdf', width = 10, height = 6)